###############################################################################
### Replication script for "International Trade, Cooperation, and Conflict: ###
### The Role of Institutions and Capabilities #################################
###############################################################################

# Analysis conducted using R version 4.0.3 (2020-10-10) on OS 10.15.6

library(stargazer)
library(countrycode)
library(ggplot2)
library(ggpubr)


# Note: this working directory must be changed to identify the "Replication" folder
setwd("FPA Replication")

final.pooled <- get(load("final.pooled.RData"))
final <- get(load("final.RData"))

# Models presented in Table 1:
### Using PR/SMD ###
model.cf.pooled <- lm(random.slope.conf.init.pooled ~ comb.reg.type + log.cinc + unitarism + v2dlencmps +
                        count.contig, data = final.pooled, na.action = na.omit)

summary(model.cf.pooled)

model.co.pooled <- lm(random.slope.coop.init.pooled ~ comb.reg.type + log.cinc + unitarism + v2dlencmps +
                        count.contig, data = final.pooled, na.action = na.omit)

summary(model.co.pooled)

### Using IDI ###
model.cf.pooled.i <- lm(random.slope.conf.init.pooled ~ IDI.IRT.mean + log.cinc + unitarism + v2dlencmps + 
                          count.contig, data = final.pooled, na.action = na.omit)

summary(model.cf.pooled.i)

model.co.pooled.i <- lm(random.slope.coop.init.pooled ~ IDI.IRT.mean + log.cinc + unitarism + v2dlencmps +
                          count.contig, data = final.pooled, na.action = na.omit)

summary(model.co.pooled.i)

# models presented in Table 2:
# yearly with clientelism

model.cf <- lm(random.slope.conf.init ~ comb.reg.type + log.cinc + unitarism + v2dlencmps + rivalry + 
                 wdi.gdpgr + count.contig + as.factor(year), data = final, na.action = na.omit)

summary(model.cf)

model.co <- lm(random.slope.coop.init ~ comb.reg.type + log.cinc + unitarism + v2dlencmps + rivalry + 
                 wdi.gdpgr + count.contig + as.factor(year), data = final, na.action = na.omit)

summary(model.co)

model.cf.i <- lm(random.slope.conf.init ~ IDI.IRT.mean + log.cinc + unitarism + v2dlencmps + rivalry + 
                   wdi.gdpgr + count.contig + as.factor(year), data = final, na.action = na.omit)

summary(model.cf.i)

model.co.i <- lm(random.slope.coop.init ~ IDI.IRT.mean + log.cinc + unitarism + v2dlencmps + rivalry + 
                   wdi.gdpgr + count.contig + as.factor(year), data = final, na.action = na.omit)

summary(model.co.i)


#Figure 1:
############################
### Random slope figures ###
############################
### For pooled models ######
############################
final.pooled2 <- subset(final.pooled, year == 1993)
final.pooled2$count.1 <- rank(final.pooled2$random.slope.conf.init.pooled)
final.pooled2$country.name <- countrycode(final.pooled2$ccode, "cown", "country.name")
final.pooled2 <- final.pooled2[order(final.pooled2$random.slope.conf.init.pooled),]

s1 <- ggplot(data = final.pooled2) +
  geom_pointrange(aes(x = as.factor(count.1), y = random.slope.conf.init.pooled, 
                      ymin = random.slope.conf.init.pooled - 1.96 * random.slope.conf.init.se.pooled, 
                      ymax = random.slope.conf.init.pooled + 1.96 * random.slope.conf.init.se.pooled),
                  color = "dark red") +
  geom_hline(aes(yintercept = 0)) +
  coord_flip() +
  theme_bw() +
  ylab("Random trade/gdp slope - conflict initiaton") +
  xlab("Country") +
  theme(axis.text.y = element_text(hjust=0, size=rel(0.65)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank()) + 
  scale_x_discrete(labels = final.pooled2$country.name)

final.pooled2 <- subset(final.pooled, year == 1993)
final.pooled2$count.1 <- rank(final.pooled2$random.slope.coop.init.pooled)
final.pooled2$country.name <- countrycode(final.pooled2$ccode, "cown", "country.name")
final.pooled2 <- final.pooled2[order(final.pooled2$random.slope.coop.init.pooled),]

s2 <- ggplot(data = final.pooled2) +
  geom_pointrange(aes(x = as.factor(count.1), y = random.slope.coop.init.pooled, 
                      ymin = random.slope.coop.init.pooled - 1.96 * random.slope.coop.init.se.pooled, 
                      ymax = random.slope.coop.init.pooled + 1.96 * random.slope.coop.init.se.pooled),
                  color = "blue") +
  geom_hline(aes(yintercept = 0)) +
  coord_flip() +
  theme_bw() +
  ylab("Random trade/gdp slope - cooperation initiaton") +
  xlab("Country") +
  theme(axis.text.y = element_text(hjust=0, size=rel(0.65)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank()) + 
  scale_x_discrete(labels = final.pooled2$country.name)

ggarrange(s1, s2)

ggsave(file="Figure1.pdf", width=10, height=8)


#Figure 2
###################################################
### predicted slope figures for major countries ###
###################################################

for.yearly.plots.setup <- final

for.yearly.plots.setup$country.name <- countrycode(for.yearly.plots.setup$ccode, "cown", "country.name")

for.yearly.plots2 <- subset(for.yearly.plots.setup, country.name == "United States" | country.name == "India" | country.name == "Russia" | country.name == "United Kingdom" | country.name == "Japan" | 
                              country.name == "Australia" | country.name == "Finland" | country.name=="Burkina Faso" | country.name=="Cameroon" | country.name=="Norway" | 
                              country.name=="Canada" | country.name == "Switzerland" | country.name == "Luxembourg" | country.name == "Belgium" | country.name == "Portugal" | 
                              country.name == "Denmark" | country.name == "Costa Rica" | country.name == "Belgium" | country.name == "South Korea" | country.name ==  "Iceland" | 
                              country.name == "Finland" | country.name == "New Zealand" | country.name == "South Korea" | country.name == "Italy" | country.name == "Greece" | 
                              country.name == "South Africa" | country.name == "Cuba" | country.name == "Columbia" | country.name == "Spain",
                            select = c("ccode", "country.name", "year", "random.slope.conf.init", "random.slope.conf.init.se", "random.slope.coop.init", "random.slope.coop.init.se"))

ggplot(for.yearly.plots2, aes(x=year)) + 
  geom_errorbar(aes(ymin=random.slope.conf.init - 1.96 * random.slope.conf.init.se
                    ,ymax=random.slope.conf.init + 1.96 * random.slope.conf.init.se), 
                width=0.1, color = "dark red") +
  geom_point(aes(y=random.slope.conf.init), size=1, shape=21, fill = "dark red") +
  geom_errorbar(aes(ymin=random.slope.coop.init - 1.96 * random.slope.coop.init.se
                    ,ymax=random.slope.coop.init + 1.96 * random.slope.coop.init.se), 
                width=0.1, color = "blue") +
  geom_point(aes(y=random.slope.coop.init), size=1, shape=21, fill = "blue") +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() + ylab("Estimated slope") +
  facet_wrap(~country.name, ncol=5) + 
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position="bottom")

ggsave(file="Figure2.pdf", width=10, height=8)


#####################################
### Predicted slopes for analysis ###
#####################################

#Figure 3:
### 3-category pooled ###
for.pred.medcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                  unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  v2dlencmps = c(median(v2dlencmps[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  log.cinc = quantile(log.cinc, probs = 0.5, na.rm = TRUE, names = FALSE),
                                                  count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE))
))

pred.cf2 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                       predict(model.cf.pooled, for.pred.medcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 50th pct")

pred.co2 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                       predict(model.co.pooled, for.pred.medcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 50th pct")

for.pred.lowcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                  unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  v2dlencmps = c(median(v2dlencmps[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  log.cinc = quantile(log.cinc, probs = 0.25, na.rm = TRUE, names = FALSE),
                                                  count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE))
))

pred.cf2.25 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.cf.pooled, for.pred.lowcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 25th pct")

pred.co2.25 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.co.pooled, for.pred.lowcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 25th pct")

for.pred.highcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                   unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   v2dlencmps = c(median(v2dlencmps[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   log.cinc = quantile(log.cinc, probs = 0.75, na.rm = TRUE, names = FALSE),
                                                   count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE))
))

pred.cf2.75 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.cf.pooled, for.pred.highcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 75th pct")

pred.co2.75 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.co.pooled, for.pred.highcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 75th pct")

pred.combined.pow <- rbind(pred.cf2.25, pred.co2.25, pred.cf2, pred.co2, pred.cf2.75, pred.co2.75)

pr.pooled <- ggplot(pred.combined.pow, aes(x=as.factor(pow), color = dv)) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=0, size=0.5, position = position_dodge(.1)) +
  geom_point(aes(y=fit), size=1, position = position_dodge(.1)) +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() + ggtitle("3-category regime type/electoral system (all states)") +
  ylab("Prediction") + xlab(NULL) + facet_wrap(~Regime) +
  scale_colour_manual(name = "", values = c("dark red", "blue")) +
  scale_x_discrete(labels = c("25th CINC", "50th CINC", "75th CINC")) +
  #                            "Maj: 25th CINC", "Maj: 50th CINC", "Maj: 75th CINC",
  #                            "Prop: 25th CINC", "Prop: 50th CINC", "Prop: 75th CINC")) +
  # guides(color=guide_legend(nrow=2,byrow=TRUE)) + 
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position="bottom")

### IDI figures - pooled ###
for.pred.cf2.i.25 <- with(final.pooled, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                                   unitarism = median(unitarism, na.rm = TRUE),
                                                   v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                                   log.cinc = quantile(log.cinc, probs = 0.25, na.rm = TRUE, names = FALSE),
                                                   count.contig = median(count.contig, na.rm = TRUE)
))

pred.cf2.i.25 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.pooled.i, for.pred.cf2.i.25, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=25th pct")

pred.co2.i.25 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.pooled.i, for.pred.cf2.i.25, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=25th pct")

for.pred.cf2.i.50 <- with(final.pooled, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                                   unitarism = median(unitarism, na.rm = TRUE),
                                                   v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                                   log.cinc = quantile(log.cinc, probs = 0.5, na.rm = TRUE, names = FALSE),
                                                   count.contig = median(count.contig, na.rm = TRUE)
))

pred.cf2.i.50 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.pooled.i, for.pred.cf2.i.50, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=50th pct")

pred.co2.i.50 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.pooled.i, for.pred.cf2.i.50, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=50th pct")


for.pred.cf2.i.75 <- with(final.pooled, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                                   unitarism = median(unitarism, na.rm = TRUE),
                                                   v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                                   log.cinc = quantile(log.cinc, probs = 0.75, na.rm = TRUE, names = FALSE),
                                                   count.contig = median(count.contig, na.rm = TRUE)
))

pred.cf2.i.75 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.pooled.i, for.pred.cf2.i.75, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=75th pct")

pred.co2.i.75 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.pooled.i, for.pred.cf2.i.75, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=75th pct")


for.idi.graph <- rbind(pred.cf2.i.25, pred.co2.i.25, pred.cf2.i.50, pred.co2.i.50, pred.cf2.i.75, pred.co2.i.75)

idi.pooled <- ggplot(for.idi.graph, aes(x=as.factor(pow), color = as.factor(dv))) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=0, size=0.5, position = position_dodge(.1)) +
  geom_point(aes(y=fit), size=1, position = position_dodge(.1)) +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() + ggtitle("Institutionalized popular inclusion (dems only)") +
  ylab("Prediction") + xlab(NULL) + facet_wrap(~ IDI) +
  scale_colour_manual(name = "", values = c("dark red","blue")) +
  # guides(color=guide_legend(nrow=2,byrow=TRUE)) +
  scale_x_discrete(labels = c("25th CINC", "50th CINC", "75th CINC")) +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position="bottom") 

pooled <- ggarrange(pr.pooled, idi.pooled, common.legend = TRUE, legend = "bottom", nrow = 2)
pooled

pdf(file = "Figure3.pdf", width=10, height=10)
pooled
dev.off()


#Figure 4
### 3-category yearly ###
for.pred.medcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                  v2dlencmps = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  log.cinc = quantile(log.cinc, probs = 0.5, na.rm = TRUE, names = FALSE),
                                                  rivalry = c(median(rivalry[comb.reg.type == "proportional"], na.rm = TRUE), median(rivalry[comb.reg.type == "majoritarian"], na.rm = TRUE), median(rivalry[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  wdi.gdpgr = c(median(wdi.gdpgr[comb.reg.type == "proportional"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "majoritarian"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  year = median(year, na.rm = TRUE)
))

pred.cf2 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                       predict(model.cf, for.pred.medcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 50th pct")

pred.co2 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                       predict(model.co, for.pred.medcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 50th pct")

for.pred.lowcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                  v2dlencmps = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  log.cinc = quantile(log.cinc, probs = 0.25, na.rm = TRUE, names = FALSE),
                                                  rivalry = c(median(rivalry[comb.reg.type == "proportional"], na.rm = TRUE), median(rivalry[comb.reg.type == "majoritarian"], na.rm = TRUE), median(rivalry[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  wdi.gdpgr = c(median(wdi.gdpgr[comb.reg.type == "proportional"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "majoritarian"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                  year = median(year, na.rm = TRUE)
))

pred.cf2.25 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.cf, for.pred.lowcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 25th pct")

pred.co2.25 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.co, for.pred.lowcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 25th pct")

for.pred.highcinc <- with(final.pooled, data.frame(comb.reg.type = c("proportional", "majoritarian", "authoritarian"),
                                                   v2dlencmps = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "majoritarian"], na.rm = TRUE), median(v2dlencmps[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   unitarism = c(median(unitarism[comb.reg.type == "proportional"], na.rm = TRUE), median(unitarism[comb.reg.type == "majoritarian"], na.rm = TRUE), median(unitarism[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   log.cinc = quantile(log.cinc, probs = 0.75, na.rm = TRUE, names = FALSE),
                                                   rivalry = c(median(rivalry[comb.reg.type == "proportional"], na.rm = TRUE), median(rivalry[comb.reg.type == "majoritarian"], na.rm = TRUE), median(rivalry[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   wdi.gdpgr = c(median(wdi.gdpgr[comb.reg.type == "proportional"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "majoritarian"], na.rm = TRUE), median(wdi.gdpgr[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   count.contig = c(median(count.contig[comb.reg.type == "proportional"], na.rm = TRUE), median(count.contig[comb.reg.type == "majoritarian"], na.rm = TRUE), median(count.contig[comb.reg.type == "authoritarian"], na.rm = TRUE)),
                                                   year = median(year, na.rm = TRUE)
))

pred.cf2.75 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.cf, for.pred.highcinc, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC = 75th pct")

pred.co2.75 <- data.frame(Regime = c("Proportional", "Majoritarian", "Authoritarian"), 
                          predict(model.co, for.pred.highcinc, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC = 75th pct")

pred.combined.pow <- rbind(pred.cf2.25, pred.co2.25, pred.cf2, pred.co2, pred.cf2.75, pred.co2.75)

pr.yearly <- ggplot(pred.combined.pow, aes(x=pow, color = dv)) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=0, size=0.5, position = position_dodge(.1)) +
  geom_point(aes(y=fit), size=1, position = position_dodge(.1)) +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() + ggtitle("3-category regime type/electoral system (all states)") +
  ylab("Prediction") + xlab(NULL) + facet_wrap(~ Regime) +
  scale_colour_manual(name = "", values = c("dark red", "blue")) +
  scale_x_discrete(labels = c("25th CINC", "50th CINC", "75th CINC")) +
  # guides(color=guide_legend(nrow=2,byrow=TRUE)) + 
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position="bottom")

### IDI yearly ###
for.pred.cf2.i.25 <- with(final, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                            v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                            unitarism = median(unitarism, na.rm = TRUE),
                                            log.cinc = quantile(log.cinc, probs = 0.25, na.rm = TRUE, names = FALSE),
                                            rivalry = median(rivalry, na.rm = TRUE),
                                            wdi.gdpgr = median(wdi.gdpgr, na.rm = TRUE),
                                            count.contig = median(count.contig, na.rm = TRUE),
                                            year = median(year, na.rm = TRUE)
))

pred.cf2.i.25 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.i, for.pred.cf2.i.25, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=25th pct")

pred.co2.i.25 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.i, for.pred.cf2.i.25, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=25th pct")

for.pred.cf2.i.50 <- with(final, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                            v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                            unitarism = median(unitarism, na.rm = TRUE),
                                            log.cinc = quantile(log.cinc, probs = 0.5, na.rm = TRUE, names = FALSE),
                                            rivalry = median(rivalry, na.rm = TRUE),
                                            wdi.gdpgr = median(wdi.gdpgr, na.rm = TRUE),
                                            count.contig = median(count.contig, na.rm = TRUE),
                                            year = median(year, na.rm = TRUE)
))

pred.cf2.i.50 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.i, for.pred.cf2.i.50, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=50th pct")

pred.co2.i.50 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.i, for.pred.cf2.i.50, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=50th pct")

for.pred.cf2.i.75 <- with(final, data.frame(IDI.IRT.mean = quantile(IDI.IRT.mean, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE),
                                            v2dlencmps = median(v2dlencmps, na.rm = TRUE),
                                            unitarism = median(unitarism, na.rm = TRUE),
                                            log.cinc = quantile(log.cinc, probs = 0.75, na.rm = TRUE, names = FALSE),
                                            rivalry = median(rivalry, na.rm = TRUE),
                                            wdi.gdpgr = median(wdi.gdpgr, na.rm = TRUE),
                                            count.contig = median(count.contig, na.rm = TRUE),
                                            year = median(year, na.rm = TRUE)
))

pred.cf2.i.75 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.cf.i, for.pred.cf2.i.75, type = "response", interval = "confidence"), dv = "Trade-Conflict", pow = "CINC=75th pct")

pred.co2.i.75 <- data.frame(IDI = c("IDI=5th pct", "IDI=95th pct"), 
                            predict(model.co.i, for.pred.cf2.i.75, type = "response", interval = "confidence"), dv = "Trade-Cooperation", pow = "CINC=75th pct")

for.idi.graph <- rbind(pred.cf2.i.25, pred.co2.i.25, pred.cf2.i.50, pred.co2.i.50, pred.cf2.i.75, pred.co2.i.75)

idi.yearly <- ggplot(for.idi.graph, aes(x=as.factor(pow), color = as.factor(dv))) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=0, size=0.5, position = position_dodge(.1)) +
  geom_point(aes(y=fit), size=1, position = position_dodge(.1)) +
  geom_hline(aes(yintercept = 0)) +
  theme_bw() + ggtitle("Institutionalized popular inclusion (dems only)") +
  ylab("Prediction") + xlab(NULL) + facet_wrap(~ IDI) +
  scale_colour_manual(name = "", values = c("dark red","blue")) +
  # guides(color=guide_legend(nrow=2,byrow=TRUE)) +
  scale_x_discrete(labels = c("25th CINC", "50th CINC", "75th CINC")) +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", linetype="dotted"),
        panel.grid.major.x = element_line(size=0.25, colour="grey80", linetype="dotted"),              
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position="bottom")

yearly <- ggarrange(pr.yearly, idi.yearly, common.legend = TRUE, legend = "bottom", nrow = 2)
yearly

pdf(file = "Figure4.pdf", width=10, height=10)
yearly
dev.off()

# Tables from paper
coef.table1 <- stargazer(model.cf.pooled, model.cf.pooled.i, model.co.pooled, model.co.pooled.i,
                         column.labels = c("Model 1", "Model 2", "Model 3", "Model 4"),
                         model.names=FALSE, style="aer",
                         ci=TRUE, p.auto=TRUE,
                         dep.var.labels.include=TRUE,
                         dep.var.labels = c("Trade-Conflict", "Trade-Cooperation"),
                         table.placement="htbp",
                         font.size="scriptsize",
                         model.numbers=FALSE, #float.env="sidewaystable",
                         notes.label=" ",
                         digits=2, df=FALSE, omit.stat=c("rsq", "aic"),
                         order = c(2, 1),
                         covariate.labels = c("Democracy-PR", "Authoritarian", "Institutional popular inclusion", "log CINC",
                                              "Unitary govt", "Public goods index", "Count of contiguous states"),
                         no.space=TRUE,
                         title = "Coefficients and 95 percent confidence bounds examining structural determinants of the relationship between trade on politics, pooled models.",
                         label = "table:coefficients1",
                         notes.align = "r",
                         notes = c("* p less than 0.1, ** p less than  0.05, *** p less than  0.01"),
                         star.cutoffs = c(0.1, 0.05, 0.01),
                         omit = c("year"),
                         notes.append = FALSE)

write(x=coef.table1, file="Table1.tex")


coef.table2 <- stargazer(model.cf, model.cf.i, model.co, model.co.i,
                         column.labels = c("Model 5", "Model 6", "Model 7", "Model 8"),
                         model.names=FALSE, style="aer",
                         ci=TRUE, p.auto=TRUE,
                         dep.var.labels.include=TRUE,
                         dep.var.labels = c("Trade-Conflict", "Trade-Cooperation"),
                         table.placement="htbp",
                         font.size="scriptsize",
                         model.numbers=FALSE, #float.env="sidewaystable",
                         notes.label=" ",
                         digits=2, df=FALSE, omit.stat=c("rsq", "aic"),
                         order = c(2, 1),
                         covariate.labels = c("Democracy-PR", "Authoritarian", "Institutional popular inclusion", "log CINC",
                                              "Unitary govt", "Public goods index", "Active rivalry", "GDP growth", "Count of contiguous states"),
                         no.space=TRUE,
                         title = "Coefficients and 95 percent confidence bounds examining structural determinants of the relationship between trade on politics, yearly models. Note:
                              all models include year fixed effects (coefficients not presented).",
                         label = "table:coefficients2",
                         notes.align = "r",
                         notes = c("* p less than 0.1, ** p less than  0.05, *** p less than  0.01"),
                         star.cutoffs = c(0.05, 0.01, 0.001),
                         omit = c("year"),
                         notes.append = FALSE)

write(x=coef.table2, file="Table2.tex")

